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On Numerically Accurate Finite Element Solutions in the Fully Plastic Range 

by 

+ 

J * C. Nagtegaal , D. M. Parks, and J. R. Rice 
Bivision of Engineering, Brown University, Providence, R. I., USA 

Abstract 

It is often found that tangent-stiffness finite element solutions for 
elastic-plastic materials exhibit much too stiff a response in the fully plastic 
range* This is most striking for the perfectly plastic material idealization, 
in which case a limit load exists within conventional small displacement gradient 
assumptions. However, finite element solutions often exceed the limit load by 
substantial amounts, and in some cases have no limit load at all. It is shown 
that a cause of this inaccuracy is that incremental deformation fields of typi- 
cal two and three-dimensional finite elements are highly constrained at or near 
limit load. This is shown to enforce unreasonable kinematic constraints on the 
modes of deformation which assemblages of elements are capable of exhibiting. 

A general criterion for testing a mesh with topologically similar repeat units 
is given, and the analysis shows that only a few conventional element types and 
arrangements are, or can be made suitable for computations in the fully plastic 
range. Further, a new variational principle, which can easily and simply be 
incorporated into an existing finite element program, is presented. This allows 
accurate commutations to be made even for element designs that would not normally 
be suitable. Numerical results are given for three plane strain problems, namely 
pure bending of a beam, a thick -walled tube under pressure, and a deep double 
edge cracked tensile specimen. These illustrate the effects of various element 
designs and of the new variational procedure. An appendix extends the discussion 
to elastic-plastic computation at finite strain. 
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1. Introduction 
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In recent years, the finite element method has been employed for analysis 
of structures exhibiting elastic-plastic material behavior [■!]* and many successful 
applications have been made (see C‘2] for a comparison of possible forms of the 
required incremental analysis). However, with the application of the method, 
at least in its tangent-stiffness form, to problems of plane strain, axisymmetric 
and three-dimensional problems, inaccurate results are often obtained. The problem 
is most clearly demonstrated for structures of ideally plastic material. Then a 
limit load exists, which can sometimes be calculated exactly or bounded from above 
by the kinematical theorem, while the finite element solution often seems to exhibit 
no limit load at all, but rather a steadily rising load-displacement curve attaining 
values far in excess of the true limit load. Similarly, for strain-hardening 
materials having typical values for the plastic tangent modulus of two or more 
orders of magnitude less than the elastic modulus, the finite element solution 
exhibits an artificially high terminal slope of the load-deflection curve in the 
fully plastic range. 

Later we shall see examples of this for a beam in pure bending and for a 

punch problem. In the case of a plane strain extrusion problem [3], when a 

- - -3 

bilinear constitutive law with slight plastic work -hardening (da/de^ = H = 4.4x10 E 
E , the Young's modulus) was incorporated, a finite element extrusion pressure of 
1.5 times a slip-line solution (using the same yield stress as the linearly har- 
dening model) was attained with overall billet displacements only of order six 
times the displacement obtained at the slip line limit load. This indicates that 
the computed stiffness in the fully plastic range far exceeds what would be ex- 
pected for the small hardening modulus used. 


In some cases, however, the correct limit load is obtained. Fig. 1 
illustrates some of the peculiarities which have been encountered by the authors 
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and co-workers in elastic-perfectly plastic analysis of plane strain bodies 
containing cracks. The finite element meshes shown in fig. 2 were used to 
solve the single-edge notched plane strain tensile specimen subject to both 
uniform end displacement £4} (i.e., no end rotation) and the static equivalent 
of a mid-ligament concentrated force. These two loading conditions have the same 
analytical limit load, namely, twice the yield stress in shear on the net section. 
As can be seen in fig. 1, the overall load-deflection curves of the two loading 
conditions are quite different in the fully plastic regime. The solution for 
the mid-ligament loading clearly indicates the correct limit load, while the 
constant end displacement solution exceeds limit load, with the load continuing 
to rise at a roughly constant rate. 

A cause of these problems relates to the fact that the deformation state 
of an elastic-perfectly plastic material is highly constrained at limit load; 
for the usual material idealization, deformation increments at limit load will 
be strictly incompressible . In the usual finite element formulation, in terms 
of kinematically admissible displacement fields, the same condition will have 
to be satisfied. 

In particular, a tangent-stiffness finite element solution satisfies the 
incremental virtual work principle 


T . u . dS 
l i 


I 

elem. 


L, 

elem. 


* o 

o . . c . . dV , 
lj lj elem 


(l.n 


precisely, where cr. . is the stress rate following from the prescribed con- 

* 3 

stitutive law in terms of the current stress a., and strain rate c. . within 

13 13 

each element. This stress rate can be written as 
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a, . = s . , + — <$.. a, , 

13 ij 3 13 kk 


( 1 . 2 ) 


where is the deviatoric stress increment. Since the plastic deformation 

1 • 

will be assumed purely deviatoric, the hydrostatic stress increment y can 

be expressed as 


1 

3 




(1.3) 


where 

k a E/3(l-2v) (1.4) 

a 

is the elastic bulk modulus and is the dilatational strain increment. 

Substitution of (1.2) and (1.3) in (1.1) then furnishes 


T.u.dS 
1 1 


elem. 'V 


[s. .e. . + k( e > 2 H dV 
i] 3,3 kk elem. 


(1.5) 


elem. 


* 

where e. . is the deviatoric strain increment. In the vicinity of the limit 

• ♦ 

load, the term s^.e^ will tend to vanish pointwise, but, as follovrs from 
plastic normality, will never be negative. Hence, from (1.5), this implies that 
the incremental finite element solution satisfies 


T.u.dS 5 l 
11 elem. 


J V 1 
elem. 



dV 


elem. 


( 1 . 6 ) 


The following result may thus be stated: 

In order for a limit load to exist for the discretized finite element 
model of an elastic-plastic problem, it is necessary that the elements 
be capable of deforming so that e, , = 0 pointwise throughout the ele- 
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ments, Otherwise (1,6) requires that the load-deflection curve be 
steadily rising - i.e., no limit load exists. 

In the next section it will be shown that, except for plane stress 
problems, the requirement = 0 severely constrains the class of de- 

formations of which typical finite element grids are capable. When they are 
forced to deform in such a constrained fashion, unreasonably high limit loads, 
or no limit loads at all, will result. On the other hand, even if no limit 
load is achieved in the finite element formulation, it is to be expected (and 
is, indeed, observed) that the slope of the load-deflection curve is reduced 
substantially from its initial elastic value at load levels near the theo- 
retical limit load. This, however, does not allow an accurate inference of 
the limit load. 

Even if the material is linearly work -hardening with a small hardening 
coefficient, the terminal slope of the load-deflection curve will not necessarily 
be accurately determined (as was argued previously and will be demonstrated in 
an example), since the material will deform in a nearly incompressible fashion 
in the fully plastic range. 

It is clear that the problems discussed here have a strong similarity to 
problems encountered in the analysis of incompressible fluids and rubber-like 
solids. The difference is that the material behaves incompressibly only as 
limit conditions are approached because the effective shear modulus then tends 
to vanish, whereas fluids and rubber-like solids behave incompressibly from the 
start because the bulk modulus is very large. The essential problem is the 
same in both cases, however, in the sense that the incompressibility require- 
ment puts too severe constraints on the possible deformation modes. 
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2. Analysis 

In the previous section, it has been observed that elements must be capable 
of deforming without change of volume pointwise if a limit load is to be obtained. 
It is therefore useful to investigate the constraints this enforces upon each 
element, and the effect of these constraints on the behavior of an assemblage 
of elements. 


Consider, for instance, the grid of 4-node rectangular isoparametric ele- 
ments shown in fig. 3. Within each element, the displacement increments are 
of the form 


(u) 



{a} + {b}x + {c}y + {d}xy 


( 2 . 1 ) 


where the vectors (a), {b}, etc., are expressed in terms of the nodal veloci- 
ties and coordinates. Now, the incompressibility constraint for plane strain has 
the form 


e + e 
xx yy 



= 0 


( 2 . 2 ) 


and this requires that {d} = { 0 } as well as that b + c = 0 , a total of 

x y 

three constraints. 

The fact that {d} = (0) means that at limit load each element has strain 
increments which are constant throughout the element. It is then evident from 

displacement continuity that all elements marked * in fig. 3 must have the 

* 

same value of e HX » whereas all elements marked + must have the same value 

of L * Moreover, since e = -e , e and s will be the same in all 

yy yy xx yy 

elements marked * or + . Since a similar argument can be set up for any row 
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• * 

and column, e and e must be the same in every element of the grid. 

Zx 

Clearly, this is a very unrealistic constraint. 

It is at least possible for this mesh to have a non-uniform shear e , 

xy 

but the situation is even worse in case the grid consists of arbitrary quadri- 
lateral isoparametric elements. Within an element the displacement increments 
are of the form 


{u} = {a} + {bin + (cH + (dine , 


(2.3) 


where n and £ are defined by 

{y} = ^ + + + * (2.4) 

(See, for instance, [6]). A lengthy but straightforward calculation furnishes 
the t ^ ree incompressibility constraints per element: 


b y - 

xy 

c B -by + 
x y y T x 

o 

H 

o 9 



b 6 - 

x y 

d 8 - b 6 + 
x y y x 

da = 
y x 

o s 


(2.5) 

d v - 

x Y y 

c 5 ** d v + 

x y y T x 

c y 6 x = 

o . 



solution to this system of equations 

in terms of three parameters 

A, B and 

b = 

X 

A6 + B3 
x y * 

b = 

y 

ce x - 

AB y , 


C e 

X 

Ay + By , 

x 'y 5 

c = 

y 

C ^x - 

% * 

(2.6) 

d 

X 

Ac$ + B 6 , 

x y 

d = 

y 

C5 x - 

AlS y • 



Substitution of (2.6) in (2.3) then makes it possible to eliminate n and £ 
from (2.3) and (2.4), which furnishes the displacement increments: 


r 



L_ 
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u = D + Ax + By 

( 2 . 7 ) 


U = E + Cx - Ay , 

y 


where the constants D and E do not depend on 
incompressibility constraint is to be satisfied, 
constant throughout the element. 


A, B and C . Hence, if the 
the strain increments will be 


This has catastrophic effects for an arbitrary grid. Consider, for instance, 
the grid in fig. and specify the displacement increments of the three nodes 
A, B and C . Because the displacement increments of nodes A and B define 
the extensional strain along AB in element I , and because the element must 
deform with a constant strain of zero dilatation, we must then specify the com- 
ponent of displacement increment of node E in the direction normal to AB . 
Similarly, because the displacement increments of nodes B and C determine 
the extensional strain along BC in element II , we must also specify the dis- 
placement component of node E in the direction normal to BC * The incompress- 
ibility constraints of elements I and II then specify the components of displace- 
ment increment of node E in two linearly independent (though not orthogonal) 
directions. Thus the displacement increment vector of E is fully determined, 
and, in fact, it corresponds to that increment which causes identical strain 
increments in elements I and II . The displacement increments of nodes D and 
F are then also determined. Continuation of the argument then furnishes, that 
the strain increment will be constant throughout the grid. It should be noted 
that this argument holds only if the element boundaries do not form a straight 
line through the body; over such a line the shear strain increment can be dis- 


continuous. 



Examples of the unreasonable constraints enforced by pointwise incompress- 
ibility upon displacement increment fields of other two and three-dimensional 
mesh configurations are given in Appendix I* 

We are now in a position to explain the results of fig. 1, bearing in mind 
a result of perfectly-plastic limit analysis: namely, that while the limit load 

is unique, the deformation field at limit load need not be unique. An acceptable 
limit field for both loading conditions, giving the correct limit load, is that 
of concentrated deformation on 45° lines from the crack tip to the surface, as 
shown in fig* 5a. An alternative limit field for the mid-ligament loading is 
shown in fig. 5b. This field consists of constant strain increments within the 
region bounded by the two 45° lines from the crack tip. The rest of the body is 
rigid, so the overall effect is that of a rotation of the two rigid regions about 
the crack tip. In fact, this was exactly the velocity field obtained at limit 
conditions in the finite element solution of the mid-ligament loading problem. 
Since the deforming regions do so uniformly, and the 45° lines coincide with 
element boundaries, the mesh of focused isoparametric quadrilaterals could 
readily accommodate the incompressibility constraints, hence leading to a very 
.accurate finite element prediction of the limit load. 

In contrast, however, the rotation produced by this field is incompatible 
with the constant end-displacement boundary condition of the other loading. The 
same constant strain pattern cannot develop, and an ever-worsening overestimate 
of the limit load results. 

Although the problem has now been analyzed sufficiently for four-noded iso- 
parametric quadrilateral elements, a more systematic approach is needed to find 
a solution to it. It is therefore useful to consider the matter of convergence 
of the solution if the mesh is refined. Refinement of the mesh will have two 
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opposing effects : 

1) It will increase the number of nodes, and since each node represents 

a specific number of degrees of freedom, it will increase the total number 
of degrees of freedom. 

2) On the other hand, it will increase the number of elements, and since 
each element has a certain number of incompressibility constraints enforced 
upon itself, it will increase the total number of constraints. 

It is clear that convergence will occur only if the total number of degrees of 
freedom increases faster than the total number of constraints. This furnishes a 
relatively simple criterion to check whether a mesh of a certain type will be 
adequate to furnish accurate limit loads and flow fields. It should be noted 
that the total number of constraints is not always equal to the number of ele- 
ments times the number of constraints per element. Sometimes elements can be 
arranged such that the constraints are no longer independent, an important 
example of which will be discussed in the next section. For the moment, however, 
these special cases will not be considered. 

It is thus important to determine the ratio of the total number of nodes 
to the total number of elements in the grid when it is refined. Consider a body, 
loaded in plane strain, the cross-section of which is subdivided into a mesh of 
elements of identical type. Define the nodal angles 6? of an element a as 
the inner angle formed at the node by the two adjacent element boundaries. For 
each type of element, the sum of these nodal angles is a given number. Consider, 
for instance, the six-noded triangle of fig. 6. Each of the nodal angles at the 
midpoints of the sides is equal to tt radians, whereas the sum of the angles 

e i* 6 2 ^ 6 3 also e 3 ual to radians. Hence the sum of all nodal angles > 

is equal to 4ir radians. Similarly simple calculations can be made for other 
types of elements. For generality, let us assume that the sum of the 
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nodal .angles of a particular element type is equal to nir radians. If the grid 
consists of p elements, then the sum of all nodal angles of all elements is 
equal to 

I I = P niT • (2.8) 

* -L 

a 1 

The sum of all nodal angles can also be calculated in a different way. 
Suppose the mesh has k nodes inside the body and & nodes on the boundary. 

The sum of the nodal angles around each interior node is equal to 2m radians, 
whereas the sum of the nodal angles at a boundary node is equal to radians, 

where y is a small angle due to the curvature of the boundary. In this way, 
the sum* is readily calculated as 

l l e“ = (2k + i + 2c - 4) IT (2.9) 

a i 


where c is the degree of connectivity of the body (c = 1 if simply connected, 
etc,). From (2.8) and (2.9) then follows the equality 


pn = 2k + JL + 2c - 4 , 


( 2 . 10 ) 


or, alternatively, 


p 2^1 

i = — ' + -T- + 


2c-4 


( 2 . 11 ) 


k n ' nk ' nk 

Now if the mesh is refined, p, k and £ will all increase; hence the last 
term will vanish. Moreover, if the mesh is refined uniformly, k (the number 
of interior nodes) will increase approximately as the square of £, (the number 
of boundary nodes). Hence 


lim 

k-*» 


<&> = 


2 

n 


( 2 . 12 ) 
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Consider what this means for a grid of isoparametric quadrilateral elements. 
The sum of the nodal angles in one element is equal to 2 tt radians, and hence 
n = 2 . Substitution in (2.12) then furnishes that in the limit the number of 
nodes will become equal to the number of elements. Since each node represents 
two degrees of freedom, and three incompressibility constraints are enforced 
upon each element, the total number of degrees of freedom is only 2/3 of the 
total number of constraints. Hence convergence will not in general occur, which 
reaffirms the previously obtained result. Similar derivations can be made for 
other planar elements as well as for axisymmetric elements, since (2.12) holds 
as well for this class. The results are displayed in Table I. The range of value 
of constraints per element from 6 to 8 given for the 8-noded quadrilateral depends 
on whether the sides of the element are linear or quadratic, respectively. 

For three-dimensional elements, however, the ratio of nodes to elements 
is not unique, since the sum of the solid angles of a polyhedron is not a con- 
stant. Therefore one has to determine this ratio for a specific arrangement of 
elements. For some common regular grids, and thus of any grid which is equi- 
valent, the results are displayed in Table II. It may be assumed that these 
results are fairly representative for most common types of arrangements. 

The results in Tables I and II indicate that none of the usual finite 
elements, except the 6-noded plane strain triangle and, to a much lesser degree, 
the 10-noded tetrahedron, is adequate to analyze (approximately) incompressible, 
and thus limit load, behavior. Before going into a more systematic approach to 
solve the problem, let us consider what can be obtained by a special arrange- 


ment of the elements. 



3 * Special Arrangements of Elements 

As has been mentioned previously, the total number of constraints is not 
necessarily equal to the product of the number of elements and the number of 
constraints per element. For special arrangements of elements, the total number 
of constraints may be less than this product. The obtained improvement, however, 
will be small, since it is necessary to combine a number of elements to eliminate 
one constraint. Thus this method will work only for elements for which the ratio 
of degrees of freedom to constraints in Tables I or II is close or equal to one. 

Consideration of Tables I and II with this in mind indicates that only a few 
types of elements could possibly be used for this purpose: for plane strain, the 

3-noded triangular element and perhaps the 8-noded quadrilateral element, and for 
three-dimensional problems , the 10-noded tetrahedron. Hone of the axisymmetric 
elements could possibly be used. As yet the only successful arrangement that has 
been discovered is the combination of four 3-noded triangles as shown in fig. 7. 
Note that the triangular element boundaries form a quadrilateral "element” and 
its diagonals. This combination has the special property that if the incompress- 
ibility constraint is satisfied in three of the elements, the constraint in the 
fourth element is automatically satisfied. This brings the ratio of degrees of 
freedom to constraints from 1 to 4/3, which makes the arrangement suitable for 
analysis near limit load. The proof of this special property is given below. 

Let the coordinates of the node i be x. , and let the node have a dis- 
placement increment u. . The area A of element I before the displacement 
increment has taken place is then equal to 


(3.1) 



where x and - denote the usual vector and scalar products respectively, and 
k is the unit vector perpendicular to the plane of the elements. Similarly one 


finds that after the displacement increment has taken place. 


• • 


« • 


A + A = — [(x +u -x -u )x (x.+u.-x -u )3 • k . 

2 2 ~2 n o ~o "1 ~1 r “0 ** 


(3.2) 


Subtraction of (3.1) from (3.2) furnishes 


A = 4 [(x 0 -x )x(u -u )-(x 1 -x )x(i -i )+(i -It )*(u,-u )] • k . (3.3) 

2 *2 o ~1 ~o ~1 "*"0 ~2 ~o ~2 ~o '-'l ~o ~ 

When one neglects the second order terms, the incompressibility constraint 
takes the form 


(x -X ) X ( U_ — U ) 
~2 v o -1 ~o 


(x -X ) X ( u*-u ) . 
~1 ~'o ~2 ^o 


(3.4) 


Similarly, for the elements II and III s one obtains 


) x (v-u) = (x -x ) X (u-u ) , 

/vj ~2 ~o ~2 v o ~ 3 " / o 


(3.5) 


(x -x ) ^ (u-u ) 

^4 ~3 ~o 


(x -X ) X (u.-u ) . 
~ o ~0 ~4 ~o 


(3.6) 


Now multiply (3.4) by I ^ 3 “b 0 i / I > furnishing 
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since - x q and x^ - x o have opposite directions. Similarly, one obtains 
for (3.6) 


. . i«o‘ • • 

-(x-x)x(u-u ) = ( X_“X )x(u.“U ) . 

2 ~o "3 — o i I ~ 3 ~o ~4 

x-x 
1 ~4 -o' 


(3.8) 


Substitution of (3.7) and (3.8) in (3.5) then furnishes 


~ o -“0 


X -X 
~o 


(x -X )x(u -U ) 

~ l ''O 


X--X I 
~2 “O 


X. -X 

/V Q 1 


( X * “ X_ ) x ( u, “U ) , 


(3.9) 


or, alternatively. 


x.-xj . x.-x 

~~ = (x -x )x(u ~u ) . 

Z O ^ 1 O I | 3 ~0 >v4 ^ o 


X„-X 
'*2 ~o 


X--X 

<~o 


(3.10) 


With the same argument as was used for (3.7) and (3.8) this yields 


• c 


(x.-x Wir-u ) 
~4 ~o "O 


(x-x)x(u-u) , 
~1 w O ^ 4 


(3.11) 


which proves that the incompressibility constraint in element IV is satisfied. 

It should be noted that this approach requires no modification of existing 
finite element programs. The groups of four elements have all the generality in 
application of the usual quadrilateral elements, and the calculations can be 
carried out without changes to an already existing program. The principal short- 
coming of the method is that it can be applied only to problems of plane strain, 
for which another alternative (the 6-noded triangular element) already exists. 
Even this alternative can be s lightly improved by arranging the 6 —node d elements 
similar to the 3— noded ones. The ratio of degrees of freedom to constraints then 
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improves from 4/3 to 16/11 . 

It should also be noted that the degrees of freedom available need not 
actually be incorporated into the finite element solution. For example the 
degrees of freedom corresponding to the interior node "0” of fig. 7 may be 
’’condensed out” of the master stiffness matrix and recovered later, (see, for 
example, [7]) 

As yet no method has been devised to arrange the 10-noded tetrahedra such 
that an acceptable ratio of degrees of freedom to constraints is obtained. 
Therefore a different, more general approach is needed to solve problems other 
than those of plane strain. 
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4. A Modified Variational Principle 

It has been discussed before that an absolute requirement for the existence 
of a limit load is that the dilatation increment vanishes pointwise. It 

has also been demonstrated that only a few of the conventional element types are, 
or can be made, useful to analyze problems under this constraint. Therefore it 
is desired that different elements be constructed, for which fewer constraints 
per element are sufficient to satisfy the incompressibility requirement. This 
goal can be obtained by taking care that the dilatation is governed by fewer pa- 
rameters than in the conventional elements. Clearly, this procedure may be applied 
only to higher order elements for which the conventional number of constraints per 
element exceeds unity. Theoretically it is possible to construct interpolation 
functions that have this property; in fact, if the arrangement in fig, 7 is taken 
as a single element, and the displacement of node 0 is chosen such that the di- 
latation in all subelements I through IV is the same, one has obtained such a 
function. This special field must, however, be considered as a "lucky guess", and 
a general approach to construct similarly simple fields for other elements does 
not seem to exist. 

An alternative way to solve the problem is to create a variational 
principle in which the dilatation al strain increment and the displacement incre- 
ments are present as independent variables. We then have complete freedom to 
characterize the dilatation by as many (or as few) parameters as we choose. Such 
a variational principle is akin to the Hellinger-Reissner principle [8], which 
admits the displacements and the stresses as independent variables. The validity 
of the present principle will, however, be proved here independently of the H-R 
principle. 
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Cobs ider the functional 


I = ICu,}3 = J CM’<4..) + - i J 2 )] dV - I T.u.dS T , 

J s T 


( 4 . 1 ) 


where 


e.» £ (u. .+u. .) - t 6. .vl , is the deviatoric strain increment; 

ID 2 i,3 3* 1 3 i] K,k ’ 

W f (e..) is the deviatoric rate potential [ W’ = 4 s . .e . . : 6W f = s..6e.. ] 

X D r 2 lj 13 13 

« 

$ is the (independent) dilatational strain increment; 
k is the (instantaneous) bulk modulus. 

This functional is well defined if a finite bulk modulus k exists. 

That this functional corresponds to a valid variational principle is readily 
proved in the usual way. By taking the first variations in u^ and <J> one finds 

61 = j^VV^Vk’d 7 ' J *i«V*T + / K( \,k' i)4idV * (4>2) 

b T 

and since 


o. . s s . . + 6 . . , 

ID ID ID 


this can be written 


( 4 . 3 ) 


61 


= l v V“i,j dV ' | s VV S T + / v K( \,k- i)5 + dV • 


( 4 . 4 ) 


The first two terms express the virtual work principle, and hence imply that 
a ij satisf y the continuing equilibrium equations in V and force rate boundary 
conditions on S T , in the usual manner. The last term provides the additional 
result that } = u^ ^ . 
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It is appropriate to compare the present principle with the principle 
created by Herrmann [5] for the analysis of incompressible, linear elastic materi- 
als. The functional used by Herrmann would have the incremental form 


H - H[u,h] = j M[e...Ew + 2vhe kk - v(l-2v)h 2 ]dV - ( T.u. 

* y j j j « ii 


i T * 


(4.5) 


where u is the shear modulus, v is Poisson’s ratio, and where the ’’mean 
pressure function; is given by 

♦ 

tl 

h = 2u(l+v) * (4.6) 


Apart from its specialization to linear elasticity, the fact that e , enters 

KK 

Herrmann's principle quadratically is an essential difference. Indeed, this does 
not allow for it the same simple procedures that can be based on the present 
principle, and that allow the easy adaption of existing, conventional stiffness 
programs to it. 


Before going to the application of the principle, let us determine how many 
constraints per element would be desired , in order to get an as close as possible 
approximation to continuum behavior. In the continuum, each material point repre' 
sents two (in problems of plane strain or axisymmetric problems) or three (in 
three-dimensional problems) degrees of freedom. For incompressibility, one con- 
straint is valid at each material point. Hence the ratio of degrees of freedom 
to constraints in the continuum is equal to two or three, depending on the type 
of problem. 


It seems plausible that a finite element solution to an elastic-plastic prob- 
lem would give the best overall approximation if the ratio degrees of freedom to 
constraints would be the same as in the continuum. Hence a "two" would be desired 
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in the last column of Table X and a "three" in the last column of Table II. This 
readily leads to a desired number of constraints per element. For instance * for 
the 4-noded quadrilateral elements in Table I , the desired number would 
be one* and for the 8-noded quadrilateral elements, it would be three. Similarly, 
for the 8-noded cube in Table II, one constraint would be desired. Similar desired 
numbers can be obtained for the other elements in the tables. 

Now let us, for the moment, restrict the discussion to those elements for 
which the desired number of constraints is equal to one. As has been discussed 
before, one then needs a dilatation increment governed by a single parameter. 

It is, on the other hand, necessary for convergence that at least a constant 
dilatation al strain increment be obtainable within an element (see, for instance, 
[7]). Hence, it is necessary to choose in an element a 




a 


constant . 


(4,7) 


Substitution of (4.7) in (4.1) then furnishes 


I a 


t , 

CX=1 i 


cw-a..) + 4 a ^ ik 


1 4 2 ] dV - 

2 a a 


T . u . dS . 
i i T 


(4.8) 


Variation of then yields the relation 



(4.9) 


In practically all cases, the bulk modulus k will be constant within an element; 
in problems of non-dilatational plasticity, for instance, k = k ^ is a constant, 
and may be brought outside the integral in (4.9) so that 
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a 


(4.10) 


Substitution of (4.10) in (4.8) then furnishes the modified functional 


1 = I L [" ,<; ij ) + * k k 1 i - L ¥i ds i • (4 - u) 

Cl ** JL V » ^ rp 


This functional can be simplified by defining the modified strain increment 


e . . = e, . + $6 . . 

J-D ID ID 




(4.12) 


so that the functional (4.11) can be written as 


= J f W(t..)dV - ( T.u.dS- , 
Jtl i V ID a J s ^ i i T » 


(4.13) 


where W represents the total rate potential (i.e., 6W(s) = a^^.6e_.^) , but 

evaluated for the modified strain increment e.. . It should be observed that the 

ID 

functional (4.13) may be expressed purely in terms of nodal displacement increments; 
hence no additional variables (i.e. , mean pressure) have to be used. It should also 
be noted that the functional is identical to the usual rate potential functional, 
with the exception of the relation between strain and displacement increments. As 
was mentioned before, this makes it extremely simple to adapt an already existing 
program for this method. One need only rewrite the matrix -converting nodal param- 
eters to strains in accord with (4.12), An equally simple formulation cannot be 
obtained with Herrmann's variational principle. 


It is interesting to note that the frequently employed procedure of replacing 
the radially- varying strains in axisymmetrie triangular elements by their values 
at the element centroid could be viewed as being based on a variational principle 
similar to the present one, in both cases implemented by changing the strain- 
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displacement relationship. 

For element types for which the desired numer of constraints is larger than 
one, the situation is only slightly more complicated. Consider for instance the 
8-noded quadrilateral plane strain element. For this element type, the desired 
number of constraints is equal to three. A possible choice for ^ is then 

5> = + (x-x ) 4> X + (y-y ) 4? , (4.14) 

a a o a 'o a 

* 0 * X *v 

where $ , <j> and <{r are constants, and x , y is the centroid of the 

a a a o o 

element. Following a procedure similar to the one followed before, one readily 
derives the relations 

*a " V 'a I v \,k dV a » 

CL 

<J> X = I^f (x-x)u, , dV 

a x a J v o a 

a 

= I ^ [ (y-y ) u, . dV 

a y a J v ° “ 

a 

where I and I are the moments of ine 
<x y o 

A similar procedure can be followed for the 8-noded quadrilateral axisymmetric 
element and the 20-noded cube. In each case e. . is defined analogously to (4.12) 
in terms of nodal displacements , and the incremental stiffness equations are 
derived in the usual way from (4.13). 

Finally it should be noted that only the quadrilateral or cubic elements 
have a desired number of constraints which corresponds to a simple and symmetric 
choice of the dilatational strain increments. For the triangular or tetrahedral 
elements the desired number of constraints does not make such a simple choice 
possible. If we did not already know that the six-noded plane strain triangle 


(4. IS) 


ctia around the x q and axes. 
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was suitable for incompressible analysis, we would conclude, for example, that 
the desired number of constraints is two, which does not correspond to a complete, 
symmetric development of the dilatational strain increments. 

Because the most significant errors in computed overall load- deflect ion curves 
occur in the fully plastic range, where the deformation gradients may no longer be 
infinitesimal, a rigorous finite strain finite element formulation may in some cases 
be required. An implementation of the variational principle presented here for such 
a formulation is given in Appendix XX. It may also be noted that the present varia- 
tional principle, eq. (4.1), may be put in a mean stress form by letting ' e q , 
say, be the variable, and this is suitable for strictly incompressible materials 
(*-*“») . 
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5. Numerical Examples 

In order to examine the results of the analyses of the preceding sections 
in specific numerical solutions, plane strain beam bending, the thick -walled plane 
strain tube under internal pressure, and the tensile analogue of the Prandtl punch 
problem were solved. 

For the first type of problem, the kinematic requirement of pure bending that 
plane sections remain plane was enforced, so it was necessary to examine only one 
row of elements across the beam thickness. Ten square 4-noded isoparametric ele- 
ments were placed across the thickness. Displacement boundary conditions were 
imposed, and the total moment was computed from the calculated nodal tractions. 

The material was modeled as elastic-perfectly plastic with a Poisson f s ratio of 
0.3. The problem was also solved by applying identical boundary conditions to 
a mesh of the constant dilatation isoparame tries developed in Section 4, and to. 
a mesh composed of 10 squares, each of which was composed of 4 diagonally crossed 
triangles, as presented in Section 3. 

The results are shown in fig. 8. The limit moment for plane strain, assuming 
a Mises yield condition, is simply 

2 

M limit ^ a o (2) (5.1) 

where c q is the tensile yield stress and h is the beam thickness. The beam 
curvature x » is normalized with respect to x Q » the curvature necessary to 
bring the outer beam fibers to yield. 

As can be seen, the ordinary isoparametric element does not find a limit load 
at all; all elements across the “thickness reach yield, and the moment- curvature 
curve continues to rise. The steepest curve shown in fig. 8 is the case for which 
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one stress state per element, calculated at the element centroid, is used to 
represent the stress state throughout the element. The next steepest curve 
results from calculating and storing the stresses at two thickness-direction 
Gaussian integration stations per element. 

The third steepest curve is the result of using the constant dilatation 
quadrilateral element with one stress state per element. Finally, the two 
essentially coincident curves which actually find the correct limit load are 
the results of the crossed triangles and the constant dilatation quadrilaterals, 
calculating and storing the stress state within each element at each Gaussian 
integration station. 

We may interpret these results by recalling that point-wise incompressibility 
is a necessary but not sufficient condition for a limit load to exist in the 
finite element model. In fact, it is also necessary that, at limit load, all 
deviatoric strain increments be normal (pointwise) to the yield surface. How- 
ever, it is not, in general, possible for the constant dilatation quadrilateral 
to deform in a manner so that the deviatoric strain increment, which is permitted 
to vary within the element, be pointwise normal to a single stress state on the 
yield surface. When the stress state was calculated and stored at the Gaussian 
integration stations, however, the computed deviatoric stress increments tended to 
vanish in the fully plastic range for both the regular and constant dilatation iso- 
parametrics. In fact, the difference between the terminal slopes of the moment- 
curvature relations for one vs. four stress measures per element was identical for 
the regular and constant dilatation isoparametrics. 

Thus, for this particular problem, the relative magnitudes of the errors in 
the terminal load— deflection curve associated with failure to meet the necessary 
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conditions of pointwise incompressibility and pointwise normality of deviatonc 
strain increments in isoparametric quadrilaterals are given, respectively, by the 
terminal slopes of the second and third steepest curves of fig. 8. 

The same problem was also solved for the case of linearly hardening material 
behavior, with a hardening modulus of one-hundredth the elastic modulus. In this 
problem, a constant terminal slope to the moment-curvature relation is readily 
calculable. The results are shown in fig. 9, where it is seen that the ordinary 
isoparametric element formulation exhibits much too stiff a response in the fully 
plastic range, a result discussed previously. The essentially coincident terminal 
slopes of the moment -curvature curves exhibited by both the constant dilatation 
quadrilaterals and the crossed triangles were less than one percent in error from 
the analytical terminal slope. 

The beam-bending problems were also solved with ten irregularly— shaped con- 
stant dilatation quadrilaterals across the thickness, and the results were es- 
sentially unaffected. 

The thick-walled plain strain tube under internal pressure proved a less 
dramatic example of the analysis of the preceding sections. A 10-degree sector 
of a tube with outer diameter to inner diameter ratio of two was modeled with 
five quadrilateral elements, the nodes having equal radial spacing. The problem 
was also solved with the same quadrilaterals formed from four triangles with the 
additional node lying at the intersection of the diagonals of the quadrilateral. 

i 

j • 

Radial displacement increments were prescribed at the two nodes on the inner dia- 
meter, and all boundary nodal displacements were constrained to be radial in 
direction. Pressures were inferred in the virtual work sense from the nodal forces 
at the two inner-diameter nodes. 

The trend in load-displacement curves obtained was the same as in the beam- 



-26- 


bending problem. However, the quantitative differences between loads at a 
given deformation were greatly reduced. In fact, the maximum difference in any 
computed load at any point in the deformation was less than four percent. The 
load obtained when '’steady state” conditions were reached was within a few percent 
of the known limit load, and the pressure-expansion curves agreed closely with 
that obtained by Hodge and White [9] using finite differences and presented 
graphically in [10 J. 

Finally, the problem of the plane strain deep, double-edge-notched (DEN) 
tensile specimen, shown in fig. 10, was solved using a very fine mesh of first 
regular then constant dilatation isoparame tries * This problem is the tensile 
analogue of Prandtl’s punch problem (see, for example, [103). In this case, the 
limit load in terms of the net stress on the ligament is given by 

°lim = (Z+ ’'>V /5 " 2 ' 97 °o (5 - 2) 

for the von Mises yield criterion, where <j q is again the uniaxial tension 
yield strength. 

By symmetry, only one-fourth of the DEN specimen was analyzed. The specimen 
width to ligament width ratio was W/b = 10 ♦ The loading was accomplished by 
imposing increments of constant end displacement at the upper end of the specimen* 
The load was computed both from the nodal tractions along the top row of nodes 
and from the average stresses. in the top row of elements. Both methods gave 
essentially the same results. 

Figure 10 shows the load-displacement curves as computed from the two 
finite element formulations. As can be seen, the ordinary isoparametric elements 
fail to find the correct limit load, and, in fact, the load-deflection curve con- 
tinues to rise at roughly constant rate, exceeding the limit value by approximately 
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25-6 at the point for which computation was stopped* Indeed, this large error 
has accumulated at an end displacement of only about 5 times the displacement 
given by extrapolating the linear elastic loading line to the limit load. 

The load-deflection curve obtained from the mesh of constant dilatation 
quadrilateral elements is indistinguishable from that of the regular isopara- 
metric elements up to roughly 60% of the limit load, at which point it begins 
to decrease in slope more rapidly. As the deformation proceeds into the fully 
plastic region the difference between the two load-deflection curves increases 
rapidly. At the point at which computation was stopped, the load was approxi- 
mately 3-5 less than the analytical limit load, and the tangent modulus load- 
deflection slope was roughly 0.1% of its elastic value, versus an apparently 
not decreasing 5-s for the regular isoparametric elements. 
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Table I. The effect of mesh refinement on the ratio of total degrees 
of freedom to total number of incompressibility constraints for some 
common two-dimensional finite elements. 
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Element Type and 
Arrangement 

Constraints 

Ratio 

Nodes 

Ratio 

Deg. Freedom 

Element 

Elements 

Constraints 

5 constant strain 
tetrahedra in cube; 
/i cubes in regular 

<7 - — lattice 

1 

1/5 

3/5 



8-node 

j isoparametric 

* cubes in 

-^regular lattice 

7 

. 

1 

3/7 

4 

l 5 linear strain 

tetrahedra in 
> cube ; cubes in 

" regular lattice 

4 

7/5 

21/20 

& 

o 

< 

— 20-node 
^ 6 °^^?i so parametric 
* -X- ©cubes in 

"Q- regular 
lattice 

5 16 

4 

$ 3/4 


Table II. The effect of mesh refinement on the ratio of total 
degrees of freedom to total number of incompressibility constraints 
for some common arrangements of three-dimensional finite elements. 
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Figure Titles 

1) Load- displacement curves obtained from finite element solutions of plane 
strain single edge cracked tensile specimens subjected to mid-ligament 
loading and constant end displacement boundary conditions. 

2) Finite element meshes used to solve the mid- ligament loading and constant 
end displacement single edge cracked tensile specimens. 

3) Propagation of incompressibility constraints in a mesh of rectangular four-node 
isoparametric plane strain finite elements. 

4) Mesh of arbitrarily skewed four-node isoparametric plane strain finite ele- 
ments . 

5a) Acceptable limit mechanism for plane strain single edge notched tensile speci- 
men subject to either mid- ligament loading or constant end displacement. 

b) Alternative limit mechanism for mid-ligament loading, 

6) Nodal angles for a linear strain triangular finite element. 

7) Four constant strain triangular finite elements arranged to form a quadrilat- 
eral and its diagonals . 

8) Moment- curvature curves obtained from finite element solutions of perfectly 
plastic plane strain beam bending. 
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9) Moment- curvature curves obtained from finite element solutions of linearly 
hardening plane strain beam bending. 

10) Load-displacement curves obtained from finite element solutions of a deep 
double edge notched tensile specimen. 

11) Propagation of incompressibility constraints in a rectangular mesh of singly 
skewed constant strain triangular finite elements. 

12) Three-dimensional block of rectangular eight-node isoparametric brick finite 
elements . 


1 



Figure 1. 
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Appendix I : Incompressibility Effects in Finite Element Meshes 

The condition e, , « 0 , pointwise, has been shown in the text to enforce 
KK 

the unrealistic constraint on rectangular plain strain 4-node isoparametric ele- 
ments that e and z must be the same for every element of the grid, and, 
xx yy 

if the 4-node i s op arame tries are arbitrarily skewed quadrilaterials , then the 
shear strain increment rate z must also be the same in every element* 

xy 

In this appendix s the consequences of = 0 pointwise will be examined 

for some other typical mesh configurations in both two and three dimensions to 
demonstrate the unrealistic constraints which are enforced upon admissible dis- 
placement rate fields* 

Consider the array of triangular constant strain elements shown in fig* 11 9 

• o 

and generated by single skewing of a rectangular grid. Since e + e = 0 , 

xx yy 

one can go from element to element along the band of elements marked * and 

• © 

show that z and c are the same in every element of such a band* 
xx yy 

Consider a rectangular array of 8-noded 9 three dimensional isoparametric 
elements, shown in fig. 12* The displacement rate field within such an element 
can be expressed as 


u 


u 


' V 


{a} + {b} x + {c} y + {d} z + {e> xy + {f} yz 

+ {g} zx + {h} xyz . (AI.l) 


The incompressibility constraint, namely 


• « 

e + e + 
xx yy 


c 

zz 



3y 



0 


(AI.2) 


T 
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requires that {h} - 0 and that the other terms be constrained such that the 
direct strains vary linearly within the element according to 

e *x = p ‘ + yz 

« 

£ yy s q - ax - yz (ai.3) 

* 

e zz = " + °x + 6y , 


where P > q s a , 3 , and y are constants within an element and are expressible 
in terms of nodal coordinates and displacement increment rates of a kind compat- 
ible with the incompressibility constraint * 


Now consider any inter-element interface of the kind z = constant. The 
» , * * 

angential strains and must be continuous across the interface and.* 

because = '^ e xx +£ yy^ always , so must be continuous across the inter- 
face. But from (AI.3), since does not vary with z within an element, and 

« 

since is continuous across element boundaries, we therefore have the con- 

straint that has the same value at all points along a line extending in the 

z-direction. That is, is independent of z , and, by similar reasoning, 

E xx ^ n ^ e P en ^ en "*- °f x and e is independent of y ♦ 


In fact, by further considerations of continuity, it may be shown that the 
normal strain rates throughout the entire array are given by equations of the 
kind 

^xx " B(y) ' C(z) 

£yy = C(z) - A(x) (AI.4) 

G zz = “ B W » 

where A(x) , B(y) , and C(z) are each continuous functions of their respective 
arguments, varying linearly within elements. 
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To see how severely this constrains admissible incremental deformation 

fields, suppose that the plane z = 0 corresponds to a fixed boundary on which 

all displacements must be zero. Then t and c vanish on this plane and, 

xx yy r ’ 

by the argument advanced earlier, so does e . Hence, from (AI.4), 

z z 

A(x) - B(y) = 0 

c(0) - A(x) a 0 (AI.5) 

B(y) - C( 0) = 0 , 

which requires that A(x) = B(y) = C(0) , and hence the most general admissible 
deformation field for the entire array of elements is 


e 

yy 

• 

e 

zz 



0 . 


C(z) - C(0) 


(AI.6) 


If we further suppose that the plane x = 0 also coincides with a fixed 
boundary, then 


e 


xx 


£ 

yy 


e 


zz 


0 


(AI.7) 


throughout * 


f 
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Appendix II: Generalization to Finite Deformations 

For the class of' elastic-plastic materials deformed at finite strain, which 
admit the existence of a constitutive rate potential, Hill [11] has shown that 
the solution of boundary value problems can be obtained from the following varia- 
tional principle, and an Eulerian finite-element formulation for problems of 
large plastic flow has recently been based upon it [12,13]: 


[ U(D )dV - i ( a |~2D D - f T.v.dSL 

I v 2 J y ij L ik k j 8>x^ J Jg i 


- 0 . 

(AII t l) 


Here V is the current (deformed) volume, and S is the portion of the current 
boundary on which nominal traction rates 1\ , based on the current state as the 
reference state, are prescribed. The rate of deformation tensor is D.. , a.. 


9v. 




13 


w » i 

is the Cauchy (true) stress tensor, and is the velocity gradient with respect 

j 

to current coordinates. The function U(D,.) is homogeneous of degree 2 in D. . 

13 i j 

and has the property that 


V _ 3U v* 

t . . - 1 r ~ ~ D « 

13 3D.. ^ijk£ k£ 

1 J 


(All. 2) 


Here are the instantaneous moduli appropriate to the adopted measure of 

stress rate, and are possibly dependent on the direction of D. , . The stress rate 

13 

. V 

itself, t.. , is the Jaumann, or co-rot ational, rate of Kirchhoff stress, based on 
a reference state that has been chosen to coincide instantaneously with the current 
state, and Kirchhoff stress t.. is defined as 0 .. times the ratio of reference 

V 

state to current material density, so that t,. = a., instantaneously . Note that 

13 1 3 

the existence of U requires that moduli have the symmetry 



and this ensures a symmetric incremental stiffness matrix for the corresponding 
finite element equations. The class of materials which admit rate potentials 
include all elastic materials for which a strain energy function exists, and also 
all elastic-plastic materials which satisfy the normality rule when phrased in 
terms of work conjugate stress and finite strain measures [14]. 

A generalization of the Prandtl-Reuss equations for finite deformation [12,13] 
which admits the rate potential U is obtained by relating the Jaumann rate of 
Kirchhoff stress to the rate of deformation tensor in the usual Prandtl-Reuss 
form: 


D? . 
12 


D?. 

13 


1+v V v , V 
- r r— t . . - — o . . t, 

E 12 E lj kk 

V 

9 T lj T )U T )c* 

4ht 2 


and 


(All. 4) 


D.. = D?. + D?. 

12 13 13 

-2 3 

for continued plastic loading. Here x = ~ and h is the slope of the 

stress- logarithmic plastic strain curve for simple tension. Note that in this 
case, and for other plastically incompressible materials, ?.. differs from 

13 

V 

only by the term (l-2v)o „o kk /E , and this is negligible by comparison for 
almost all .technological applications of the theory to metals. 

Thus equations (All. 4) can be inverted, and the resulting Jaumann rates of 


Kirchhoff stress separated into hydrostatic and deviatoric parts to give 
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• & [”« 
E 

= 1727 D kk 


- a 


>***&■ J 


s 3<D, 


kk 


(All. 5) 


where a 


1 if at yield and > 0 

0 otherwise . 


Thus the deviatoric and hydrostatic Jaumann rates of Kirchhoff stress have 
functional dependence only on the deviatoric and hydrostatic parts of the rate 
of deformation tensor, respectively. 


From these relations, the volume integral of the rate potential 
U s y in (AII.l) can be separated into deviatoric and hydrostatic parts to 

obtain 

[ i?..D..dV = I ?! .D!. + - kD* IdV . (All. 6) 

Jv 2 13 13 J„ L 2 13 13 2 kk J 


Because this integrand has the same functional form as the conventional 


small-strain rate potential for Prandtl-Reuss materials, the variational principle 

1 2 

of Section 4 can be readily generalized by replacing the term ~ icD^ in (All. 6) 

* 1 * 2 

with - as in (4.1). As in the text, independent variation of the 

+ * 

velocities and of 4> furnishes $ = , the dilatation rate. 


The finite element implementation of the finite strain version of the varia- 
tional principle of Section 4 then follows straightforwardly, with the total 
element stiffness matrix being the sum of the stiffness corresponding to the modi 
fied rate potential and the "initial stress" stiffness, which corresponds to the 
second volume integral of (All. 1). 



